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' Using a set of general methods developed by Krotov [A. I. Konnov and V. A. Krotov, Automation 

and Remote Control, 60, 1427 (1999)], we extend the capabilities of Optimal Control Theory to 
the Nonlinear Schrodinger Equation (NLSE). The paper begins with a general review of the Krotov 

Q approach to optimization. Although the linearized version of the method is sufficient for the linear 
I Schrodinger equation, the full flexibility of the general method is required for treatment of the 

1^-^ . nonlinear Schrodinger equation. Formal equations for the optimization of the NLSE, as well as a 

concrete algorithm are presented. As an illustration, we consider a Bose-Einstein condensate initially 
at rest in a harmonic trap. A phase develops across the BEC when an optical lattice potential is 
turned on. The goal is to counter this effect and keep the phase flat by adjusting the trap strength. 
' The problem is formulated in the language of Optimal Control Theory (OCT) and solved using the 

O ' above methodology. To our knowledge, this is the flrst rigorous application of OCT to the Nonlinear 

^ , Schrodinger equation, a capability that is bound to have numerous other applications. 
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I. INTRODUCTION 



■ In recent years much progress has been made in the use of Optimal Control Theory (OCT), to coherently control 
O I quantum mechanical systems governed by the Schrodinger equation. Such systems include controlled manipulations 
I ^ I . of molecular wave packets ]|[ |^ , high harmonic generation jj] and realization of quantum computing algorithms 



PsJ In this paper for the first time the application of OCT is extended in a systematic way to systems governed by 

the NLSE, such as solitons in fiber optics and Bose-Einstien condesates (EEC's) in atomic physics. We begin with a 
' general description of the Krotov iterative method |^ . We describe first its application to quantum systems governed 
by the linear Schrodinger equation and then show how a generalized version of this method Q can be used to treat 
non-linear problems. 

Finally we consider a concrete problem governed by the NLSE, namely a BEC evolving under the Gross-Pitaevskii 
fSj ' equation. The use of a BEC as a realization of quantum computing is widely being considered, as this is a macroscopic 
, entity which nevertheless behaves quantum mechanically. The fact that a BEC carries a definite phase which can 
^ ' be manipulated and controlled is a striking manifestation of this quality. For many computation applications it is 
] desirable to split the BEC up into localized pieces each of which can then be viewed as a quantum bit and manipulated 
' as such. This is achieved by the switching on of an optical lattice potential; however the switching on of the optical 
' I lattice causes a phase to accumulate across the BEC which is undesirable for use in computing applications. It is the 
cancellation of this effect which is the goal of this model problem. 

The outline of the paper is as follows: In section |^ the Krotov method is reviewed. Its application to linear problems 



in general and the Schrodinger equation in particular are discussed. Section [II deals with the application of OCT to 
the NLSE problems and demonstrates this by solving the BEC problem mentioned above. Finally in section IV we 
conclude and suggest further applications. 



II. KROTOV METHOD OF OPTIMIZATION 

A. Description of problem 

Consider a state of some system which can be defined by a vector of variables and which is controlled by a set of 
variables u, through the state equations of motion 

i; = f{t,^,u). (1) 



The initial value of t/j, "0(0) = V'o, is fixed but evolves over time to some final value 'ip{T) = ipx- The history of 
evolving state vectors is called the state trajectory and the history of control input is termed the control history or 
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just the control ^. 

Given the trajectory tp{t) and the control u{t), we define a 'process' w — {ip{t), u{t)) as a pah of histories ip{t), u{t) 
satisfying eq. (|f|). We can now define the objective functional on the process w: 

J[w] = F{^/j{T)) + r fit, m, Ht))dt, (2) 



"'0 

where F{ip{T)) and f{t,tp{t),u{t)) are general functions that represent the dependence of J on the terminal and 
intermediate time values of V' respectively. It is required to find a process w for which the objective obtains its smallest 
value. 

B. Utility constructs and definitions 

For a continuously differentiable scalar function (j){t, ip), we define the following functional 

L[w;(f,]^G{^T)- f R{t,m,u{t))dt-m^o), (3) 
Jo 

where 

G(V't) = F(V^t) + '^(T,Vt), (4) 
R{t,i:,u) = ||/(i,V;,u)-/°(i,0,w) + ^<^(t,^). (5) 

The functions R and G are designed to separate out the dependence of L[w,(j)\ on the final time and intermediate 
time respectively. 

It can be shown that for any function and process w — {'ip{t),u{t)), L[w, 0] — J[w]. The derivation goes as follows: 

L[w;4,] = G{^t)- f R{t,4,{t),u{t))dt~<j}{Q,i^o) 
Jo 

= F{^t) + HT,M- I |||/(<,V',")-/°(i,^,^) + ^|di-0(O,^o) 
= F{^t) + 0(T, £ (^^^ + ^~fit,ij,u)^dt~cp{0,i^o) 

= F{',Pt)+ f fit,^,u)dt^J[w]. (6) 
Jo 

Obviously therefore, minimizing L[w](j)\ for any minimizes J[w], and minimizing L[w,(j3\ can be achieved by sepa- 
rately minimizing G{iI>t) and maximizing R{t,'ip,u). 

It is convenient for later reference to define the function H through the following relation 

R{t,ij,u)^H{t,i^,u,^^^) + ^^cl,{t,^), (7) 

where 

Hit, ij, u,p) = pft, ^, u) - fit, ^, u). (8) 

Note the extra parameter in H denoted p, which emphasizes that ip and ^ should be treated as independent variables, 
with respect to H . 



Here and throughout the following, vector treatment of the appropriate variables is assumed although all vector notation is omitted 
to avoid congestion. Multiplication of vector variables is therefore to be understood as a dot product. An indexed notation of vector 
components will be used only when unavoidable for clarity and summation convention will then be implied. 
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C. An iterative algoritlim 



We now return to our main goal and describe the Krotov iterative method for finding a sequence of processes {ws} 
which monotonically decrease the value of the objective J[w] ||]. The central idea is that as we have complete 
freedom in choosing the potential construct such that our current estimate of the state history will 

maximize L[w; 0], and so become the worst of all possible histories. We are then free to find a new estimate for the 
control u{t) which will minimize L[w; (j)] with respect to its explicit dependence on u{t), without worrying about the 
effect of u{t) on L[w; cj)] through the change of ip{t), as that can only be improved. 

We begin by taking an arbitrary control history u^{t) and the corresponding state trajectory ■^''^{t) which constitute 
together a process 

1. We first construct a function (f>{t, ip) such that L[w; (j>] is a maximum with respect to il>{t) at the point ui". This 
is equivalent to the following two conditions: 



= mini?(<,?/>(i),u"(i)) 



max G(V't), 



(9) 
(10) 



where the functions R and G are calculated using the new (f){t, tp)- In other words we choose (j){t, ^) such that our 
current ijj^{t) will be the worst of all possible in minimizing the objective L[w; ip] — J[w] (maximizing i?, 

minimizing G). Any change in brought about by a new choice of u{t) will now only improve the minimization 
of J[w]. (see figure fl) 



R(^) ' 















FIG. 1: Sketch of various R{tp) determined by different choices of (pit, tp)- The construction of p{t, tp) is the one that minimizes 
R{tp) at ^p = ipo, i.e. is the worst R at the current ipQ. 



2. For (p{t,tp{t)) we find a control u{t) that maximizes H{t,ip,u, ^) and denote it by 

d(p 

u(t,ip) = arg max Hit, ip,u, 

u dip 

= arg max R(t,ip,u). 

u 

Note that the control u{t, ip) is still a function of ip. This freedom will be removed in the next step. 



(11) 



3. We require that u{t,ip) and ip(t) be consistent with each other through the equations of motion. The equation 
of motion (^ (with its initial conditions) together with the equation for the control u — u{t,ip) (|ll]), provide 
two equations for the two unknowns u and ip. These equations may be solved self consistently for u and ip(t), 
obtaining the new process w = (u, ip) . 

4. ft is now guaranteed that minimization of the objective has been improved so that J[w] < J[w^]; this completes 



the current iteration. The new w becomes a starting point for the next iteration, w —> 
can now be repeated to achieve further decrease in the objective. 

We proceed to prove that indeed the new J[w] < J[w^]. First we assert using eq. (||) that 

J[w°] - J[w] = L[w°] (p] - L[w; (p]. 



and operations 1-3 



(12) 



Also 
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L[w°- </)] - L[w; </)] = G(7A§,) - G(^t) + / {R{t, uit)) - Rit, V^" W, u'>it))} dt 

JQ 

= A1 + A2 + A3, (13) 

where 

Ai = G(^0)-G(V't), (14) 

A2 - / {R{t,m,m) ~ Rit^m^u^it))} dt, (15) 

Jo 

T 

{i?(t, u"(t)) - i?(t, V°(i), w°(i))} dt. (16) 

The nonnegativeness of A3 and Ai follow from conditions (^) and ( p^ respectively, and by eq. (^l|), the choice of a 
new control ensures the nonnegativeness of A2. This completes the proof. 

D. Construction of (j) to 1st order in 

In implementing the above iterative method the main difficulty lies in step 1. Here it is necessary to determine a 
function (f>{t,ip) that, by conditions and (p^, will ensure the absolute maximum and minimum of the functions 
R and G respectively on the trajectory tpo{t), i.e. to choose 4>{t,ip) to give the worst possible L[ipo;(j>]. A necessary 
condition for an extremum of R and G at = {ip^^u'^) is the existence of a stationary point there, but in order to 
make the conditions sufficient, it is necessary to add conditions of positivity and negativity on the second derivatives 
of R and G respectively. 

We leave the additional requirements on the second derivatives for a later section and restrict ourselves in this 
section solely to determining the conditions for a stationary point in R and G, which are as follows: ^ 

= |;/f(,,*»,»».,) + f ^0 (17) 



d4>T dipT dip' 



where^ 



'T 



^{t) = -^cP{t,iP"(t)). (19) 



In the derivation below the following delicate point should be noted: R{t, tp, u) is by definition a function of three variables, whereas 
Hit, ip, u,p) is a function of an additional argument p = dtp/dtp. In other words although p depends on ip with respect to R, this is not 
the case with respect to H where tp and p are to be taken as independent variables. In eq. ( |l7[ ) we wish to vary R with respect to ip 
which means H must be varied in both tp and p; dR/dip <-> {{dp/dtp)d/dp + d/dip}H . 

We stress that in the following definition, x(t) is a function solelvof t and is obtained by inserting the explicit dependence of ip^{t) into 
d<f>{t,ip) / dip. This explains the advance to the third line in eq. (O) where the total time derivative of d<p{t,ip)/dip\^o is translated to a 
simple time derivative of x(*)- 
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Eq. (0) and (|l|) can be restated as an equation of motion for the vector x and the definition of H[t,'ip^u,p)^ eq. 
can be used to rewrite the equation of motion for i/j, eq. (|l|), in the following compact form 

with boundary conditions x{T) — 5-7- (20) 



dip' 



•t 

ij = ^H{t,^y,x) 

ox 

with boundary conditions — ipo- (21) 

These equations constitute a so called Hamiltonian system with a Hamiltonian H{t, ijj, u, x), and the variables ip and 
X are said to be conjugate. Eq. ( |l9|) shows that the variable x represents the function (j){t, tp) to first order in ip. Eqs. 
( pO| , pl[ ) give the prescription for constructing x{t)- 

Creating the conjugate variable x to fulfill the above requirements enforces a stationary point in R and G with 
respect to ip. As explained above it is also necessary, in the general case, that the stationary point be an extremum and 
therefore that the second derivatives of G and R with respect to tp are negative and positive respectively. However for 
problems linear in ip it so happens that a stationary point is sufficient. This will be illustrated by a concrete example 
in the next section, after which we return to our main line of discussion completing the conditions for obtaining (j)(t, ■0) 
in the general nonlinear case. 



E. A linear problem and application to the Linear Schrodinger equation optimization 

Consider a problem where the equations of motion are linear in the state variable 

^ = a{u)^j, F{^PT) = h^T, f^f{u). (22) 

We proceed to show that it is sufficient to choose (j>{t,ip) = xlOV" to achieve monotonic increase in the objective at 
each iteration. The Hamiltonian for this problem by the definition, eq (||), is H{t,ip,u,x) — X^{'^)4' ~ f^{'^)^ so we 
get by applying eq. (|^ 

X = —a{u^)x, with boundary conditions x{T) = —b. (23) 

Using the above we find that 

Rit,^,u') = xa{u°)4'^f + ^^P = {xa{u'^) + ^)^-f 

= -f{u°), (24) 
G{iPt) = bipT + Xt^t ^ {b + Xt)^t = 0, (25) 
which are independent of ip. This implies that both Ai and A3 in eqs. (|l^ ) and (^) vanish. By maximizing A2 (eq. 



|15|)), the objective is guaranteed to decrease at each iteration. Note that in the linear case there is no need to check 
the second derivatives of G and R since, as R is linear in ip and we set §^|^o = 0, i? must be independent of ip. (see 
figure |l|) Therefore the control u can be made to maximize R without the resulting change in ip{t) having any effect 
on the objective. 

The above example encompasses the problem of optimization of a quantum mechanical wave function governed by 
the linear Schrodinger equation 

\i^)=~iH{um. (26) 

Some care is necessary, however, if the objective takes the form F{'iPt) = — {^Pt\P\''Pt) , which is not strictly linear in 
the state vector as in eq. (|2|). Another complication arises from the fact that the state vector \p is an element of a 
complex Hilbert space. We therefore work this problem out in full, and show that nevertheless the above choice of (p 
is sufficient in these problems just as in the linear example, due to the fact that the target projection operator P is 
positive definite. |l^ 

As above, we set (pltyip) = (xlV') + (V'lx) (where we have included the complex conjugate as an extra independent 
state variable) and thus get for the Hamiltonian of the problem 

H{t,^P,u,x) = ~i{x\HW +i{m^x) - 

= 2-s{x\Hium-\fiu), (27) 
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which using eq. (|20|) yields for the conjugate vector 

ix) = -iHHu'm. 

with boundary conditions \xt) — P\'^t)- (28) 
Reinserting this equation and its complex conjugate into the formulas for R and G we have 

R{t,i^y) - -*(x|i?(^i°)|V^)+*(V'|i?^(Olx)-A/°(0 

+(xl^) + {Mx) 

= -A/°(u") (29) 
G'('0t) = -('0t|-P|V't) + (xtIV't) + (V'tIxt) 

= -(V-tIPI^t) + (^§^|/'|^t) + (^t|P|V't)- (30) 

R is independent of '0 as above, which guarantees that A3 of eq. (|l6| ) vanishes. G is dependent on t/'t; however 
(denoting A-0 = tp — ipp) the positiveness of Ai = (AV'tI^'I Ai/jt) (eq |14|) is always guaranteed due to the positiveness 
of P. Alternatively, note that the second derivative of G, gj^^g^, = —P < is always negative due to the positiveness 

of P, assuring that the condition for a maximum of G'(V'o) is automatically met. ^ 

Another intricate point regarding these problems which is often missed, is the following. In many problems in 
quantum mechanics it so happens that the equations of motion are linear in the control variable u, namely the field. 
This means that strictly speaking there is no proper maximum in the Hamiltonian of the control system i?(i, "01 x) 
with respect to m, and stage 2 of the algorithm (eq. (pi])) cannot be properly fulfilled. This problem is often overcome 
by adding a penalty function, A/°(m) quadratic in u to the objective as implicitly indicated above. The physical 
interpretation of this construction is that placing a penalty on the fluence of the field, constrains the algorithm to 
search out the optimal direction of u rather than minimizing the objective by varying its magnitude. The price paid 
by this solution, however, is that the algorithm often exerts much effort into minimizing the superficial penalty part 
of the objective at the expense of the really required terminal part. 

An alternative way to overcome this problem is by noticing that the algorithm does not really require that the 
Hamiltonian be maximized by u at each iteration. All that is really required is that the Hamiltonian be increased by 



the new choice of u, which is enough to ensure that A2 (eq. (10)) be nonnegative. The penalty function can therefore 
be dropped and at each iteration u should be increased u u + A-i|f , where A"! is some macroscopic constant 
which can be chosen arbitrarily.^ This is not to be confused, despite the formal similarity, with the gradient methods 
where the correction to u must always be small such that its effect on any change in ■0 will remain in the linear regime. 



F. Construction of to 2nd order in tp 

As noted above, for an equation of motion non linear in the state variable, it is necessary to fulfill conditions on 
the second derivatives of R and G and therefore must be chosen to contain higher orders in 0. We therefore take 

1 

2' 

where Aipi = ipi — ip'^, and the functions Oij are to be determined such as to obtain the required extrema in R and G. 
The conditions supplementary to eq. ( p"7[ - p^ ), necessary for fulfilling eq. (^10), are the following system of differential 
inequalities: 



0(i, 0;) = x^^^ + -A0*ay(i)A0;, (31) 



^„d^ R{t,r,u °) 

dtp* dip j 



d^R > 0, d^R = Mi ZI J ^ V'j, (32) 



d^G < 0, d^G = A#,|!^MlAV.T, (33) 



* Another way to map the quantum control problems onto the linear example is by formulating the Schrodinger equation in the density 
matrix formalism of Liouville space like so: \p)) = iCh,\p)) with the objective F{ipj') = {{B\p)) where B is some target state. In this form 
the equivalence to the linear case becomes self evident. 

^ The choice of is not restricted by the algorithm, however some limit must exist through the physics and numerics of the problem. 
The A can also be adjusted at will during the optimization to improve convergence. The choice of the symbol A intentionally points to 
the similarity of the role this parameter plays with that of the strength of the penalty imposed above and to the similarity of the resulting 
forms for choice of the new u. The difference ends up being solely that the proposed method takes X~^dH/du to be a correction to the 
old field instead of an entirely new choice. 
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For the positivity and negativity of the quadratic forms (PR and cPG respectively it suffices to set 

Ri^^r^) = (5..n t^:^/ (34) 



di>*d^j ^ ' ' ' \6^>0 for alH = j 



dip^P^dipTj \-ar<0 for aU i = j 

where ai and Si are some nonnegative functions. Inserting the full form of R and G into the above equations yields a 
set of n{n + l)/2 equations of motion for the functions aij, where n is the dimension of the state vector. This means 
that the dimension of the system grows proportionately to the square of n, and for large n the expense of solving this 
system normally becomes too high. 

In it is proved that for certain classes of functionals conditions ([9|pO|) can be fulfilled by taking (p according to 
eq. ( pi|) with 

where a,P < and 7 > 0. Taking /3 < — ^^^g^^ always fulfills condition ( |35| ) and it can be shown that as 7, \a\ 00, 

for these classes of functionals, condition (^ ) is also fulfilled. The strategy then is to begin with a = and if the 
objective does not decrease, take increasingly larger 7, |a|, |/3| until a decrease in the objective is achieved. 

III. APPLICATION TO THE NLSE 
A. General formulation 

We wish to apply this algorithm to optimizing a quantum system governed by the nonlinear Schrodinger equation 

1^) = ^iHnlM'W = -i{H + m'^m. (37) 

where H — K -\- V is the usual linear Hamiltonian operator consisting of kinetic and potential parts and ^ is the 
coefficient of the additional non linear term. The objective is defined as for the linear case, as minimizing 

J = -(^tIPIVt) + dtf{u). (38) 

In realization of step 1 of the iterative method we choose 4i{t,ip) — (xlV") + (V'lx) + and find the 

Hamiltonian of the system to be 

= 2Q{x\HNLm-Xf{u), (39) 

just as in the linear case. However it must be remembered that here Hnl depends on tjj and ip* , and therefore using 
eq. ( pO| ) yields for the conjugate vector 

Ix) = -^i/(t,V°,"°,x) 

= ^^i{^y\x*)-^{H^ + 2fi\tPX)\x)- (40) 

Note that this equation differs in two respects from its linear counterpart (eq. First, it involves which 

does not cause special problems except that the vector j^") must be stored and used in propagating |x). The more 
cumbersome difficulty arises from the fact that the equation obtained is no longer linear in \x), but rather contains 
an extra term linear in |x*)(= (xD- Nevertheless the coupled equations 



(41) 
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where K = a — y{x) + 2^\ip^\'^ and h — ^{ip^Y , or in matrix form 

k^~i{K + V)Z (42) 

can be solved using the split operator method for separating the treatment of the kinetic and potential parts of the 
Hamiltonian: x(i + ^i) = e~*^''*/^e~*^''*e~'^'^*/^x(0- The kinetic evolution can be computed using the usual Fourier 
transform methods and the potential part can be evaluated by diagonalization to give, 



cos{D5t) - ae'-^'5* ib sm{DSt) 

a-D 



' ih* sm{D6t) cos{DSt) - ae-'^^* ' ' 



with D = ^ a? — Finally, introducing the complex conjugation operator C and the Fourier transform operator 
Z, this procedure yields the following formula for the numerical propagation step 

|x(t + 5t)) -i^e*T^**z|^ (^^L^cos{D5t) - acx^[iD8t) \ + i|- sin(D(5t)C'l Zte*^**i|x(t)). (44) 

Having obtained |x(i)) for all t, we now proceed to realize step 2 of the algorithm and according to eq. ([ll|), find 
for each point in time a field which maximizes i/(t, u, |^) — 23(x + \aA'\\)\ll{u)\'\li) — \ f\u). Or mathematically 

u{t,tlj) = argmaxH{t,ip,u, ——). (45) 

Step 3 of the algorithm is fulfilled simultaneously with step 2, by simultaneously propagating u and "0 such that 
each new u{t) is used directly in propagating iplt) ip{t + St). For a we use formula (pq ) and according to the 
algorithm described in the previous section we begin with ct = and increase the parameters 7, |a|, |/3| until achieving 
decrease of the objective. 



B. Application to a concrete problem 



We consider a ID Bose Einstein Condensate (BEG) confined by a harmonic trap and governed by the Gross- 
Pitaevskii equation 

\^) = -i{k + V + NUo\^P\^m, (46) 

where K, V are as above and NUq is the nonlinear atom-atom interaction strength, N being the number of atoms. 
The BEG is initially in the ground state of the trap potential and is therefore stationary. An optical lattice is then 
switched on, having the effect of separating the BEG wave packet into a series of localized pieces. The potential 
energy operator therefore takes the form V = Kx^ + S{t)Vo cos^{kx), where K is the trap constant, k is the laser field 
wave number, Vq is the lattice intensity and the switching on function of the field is denoted S{t). In applications to 
quantum computing, these localized wave packets are to represent quantum bits. However, due to the nonlinearity 
of the equations, the condensate develops a phase that varies from lattice site to lattice site (see figure ^), which is 
undesirable for quantum computing, since these algorithms assume that there is zero relative phase among the various 
single quantum bits. The problem is therefore to eliminate this phase profile by adjusting the trap strength during 
the switching on of the laser field. From the OGT perspective the trap constant K{t) is taken as the control and the 
objective is to minimize the variance of the phase of the wave packet, 9{x), at some final time T. The phase being a 
multivalued function poses problems; we therefore consider instead minimizing the variance of cos{9) = ^ '^j^^ such 
that the objective becomes 

J = (cos^ (61t)) - (cos(6't))^ 

= (V'|cos2(0t)|V')-(^|cos(0t)|V')' (47) 

Using the first part of eq. ( pO| ) we get the equations of motion for x as in eq. ( pO| ) and taking a derivative of J with 
respect to -0^ we get according to (|o|), the boundary conditions: 

XT = = -5R[^t] + i|^T|(cos(0T))(|r +3). (48) 

The Hamiltonian of the problem is H{t, V', K, x) ~ K{t){x + ^a{Aip)\x'^\4>) so that according to the above procedure 
we improve K at each iteration by K ^ K + X^^{x + ^(^{Atpjlx^lip) . 



Initial Wave Pacl<et {t'=0) 




Final wave function (t'=T) with no trap adjustment 

4 1 , 1 , , , , , 1 , ^5 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 



X=X/X^ 

FIG. 2: Wave packet at t=0,and t=T with no trap adjustments. 



C. Optimization Results 

Following PI we transform the NLSE to dimcnsionless units t' — t/tQ and x' — x/xq where a:o — ^tf — 20.3/im 
and to = = 75ms. The Thomas- Fermi radius xtf — \J'^IJ-tf /Tnujfj.g^p gives the size of the condensate and is 

defined through a chemical potential fiTF = ^T-ftNL determined by normalization of the wave function to unity. The 
wave function too is scaled V' ^/xtfiP and in order to keep the time scales of our ID model comparable to the 3D 
reality we adjust U CU by a factor C = r(2+I/2) ~ I l|l3- ^^^^ ^^^^ ~ 96.2/l(s, optical wavelength A — 589nm 
and Vq — 10.94i?fl for the final field intensity. All parameters were taken to resemble the experiments described in 
ipTf . Performing these transformations we end up finally with a dimcnsionless NLSE, 

= (- + K{t)x^ + S{t)VcosHkx) + C/IV'p) |V>, (49) 

where the trap constant K = cj^^^jptg, the field Intensity V = Voto/h = 2 x 250^ and the nonlinear coefficient 
U = |/iTi^to/^ — 1039, such that all space, time and energy quantities are now expressed in units of xq, to and H/to 
respectively. 

Initially the wave packet is in an eigenstate of the potential with trap constant Kq — 779 and is therefore static. 
The switching on function plotted in figure ^, turns on the optical potential at a quarter of the optimizing interval 
(r/4) and is constant from there to the final time. With no adaption of the trap constant a phase develops across 
the wave function as shown in figure ^ The optimization process decreases the objective monotonically as plotted in 
figure H and yields a strategy of increasing AK{t) = K{t) — Kq to achieve a flat phase at the final time, T = 1500/is. 
These striking results are shown in figure 0. 



IV. CONCLUSIONS 



Using a set of general methods developed by Krotov we have extended the capabilities of Optimal Control Theory to 
the Nonlinear Schrodinger Equation (NLSE). Although the linearized version of the method is sufficient for the linear 
Schrodinger equation, the full flexibility of the general method is required for a rigorous treatment of the nonlinear 
Schrodinger equation. Mention should be made of the interesting recent work of Hornung and de Vivie-Ricdlc [ p2[ , 
applying optimization techniques to molecule formations in a BEC, although the (j) function in that work included 
linear terms only. A parallel approach was pursued by Potting et al [O] who used genetic algorithms to control the 
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Objective decrease vs. Iteration 
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Iteration number 



FIG. 3: Objective decrease as function of iteration. 
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FIG. 4: Optimized trap strength evolution (bottom), and final wave packet at t=T (top). The fiat phase is strikingly apparent. 



momentum state of a BEC. The significance of this paper is twofold. First, both formal equations and a concrete 
and efficient algorithm were presented for optimizing the NLSE in cases where the nonlinear terms are significant. 
Second, the methodology was applied successfully to an interesting physical problem confronting the use of trapped 
Bose-Einstein condensates (BECs) for quantum computing, namely producing a constant final phase profile across the 
condensate after an optical lattice is turned on. Further work on understanding analytically the mechanism found by 
OCT is still underway. We believe the working equations developed here will have many more applications in systems 
governed by the NLSE, including both BECs and soliton fiber optics. 
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